# This script will demonstrate how to analyze the amphibian data collected as
# part of the LongTerm Ecological Monitoring Initiative
#
# Only one study area at time can only be analyzed with a script.
#
# This was programmed by Carl James Schwarz, Statistics and Actuarial Science, SFU
# cschwarz@stat.sfu.ca
#
# 2017-03-27 Separated egg mass by species. Imputed 0 counts when species not seen in a year
# 2017-02-28 First Edition
# Summary of calling protocol
# Define a survey transect along or through your wetland. …
# Surveyors visit the monitoring site(s) in spring and listen for calling males,
# recording the species and approximate number of each.
# Repeat surveys increase the probability that species will be detected.”
# Summary of visual protocol
# Surveyors monitor breeding site(s) during the active season (spring to fall),
# walking the shoreline of a wetland recording all species and life stages
# encountered. These include egg masses, tadpoles and adults.
# Repeat surveys increase the probability that species will be detected.”
# A separate analysis is done on EACH species present
# load libraries
library(car) # for testing for autocorrelation (2 libraries needed - see dwtest)
library(ggfortify) # for residual and other diagnostic plot
## Loading required package: ggplot2
library(ggplot2) # for plotting
library(lmtest) # for testing for autocorrelation
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(plyr) # for group processing
library(readxl) # for opening the Excel spreadsheets and reading off them
library(reshape2) # for melting and casting
library(lmerTest) # for the linear mixed modelling
## Loading required package: Matrix
## Loading required package: lme4
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(stringr) # string handling (like case conversion)
# Load some common functions
source("../CommonFiles/common.functions.R")
cat("\n\n ***** Amphibian Analysis - Single Site Analysis ***** \n\n")
##
##
## ***** Amphibian Analysis - Single Site Analysis *****
# get the data from the Excel work.books.
# we put the list of work books here, including the file type (xls or xlsx).
# You can put multiple stations here because the station information is included on the raw data
work.books.csv <- textConnection(
"file.name
General Survey Using Transects-amphibians_AliceLake2013.xls
General Survey Using Transects-amphibians_AliceLake2014.xls
General Survey Using Transects-amphibians_AliceLake2015.xls
")
work.books <- read.csv(work.books.csv, as.is=TRUE, strip.white=TRUE, header=TRUE)
cat("File names with the data \n")
## File names with the data
work.books
## file.name
## 1 General Survey Using Transects-amphibians_AliceLake2013.xls
## 2 General Survey Using Transects-amphibians_AliceLake2014.xls
## 3 General Survey Using Transects-amphibians_AliceLake2015.xls
# read each workbook and put all of the data together into one big data frame
amph.df <- plyr::ddply(work.books, "file.name", function(x){
file.name <- file.path("Data",x$file.name)
data <- readxl::read_excel(file.name, sheet="General Survey")
data
})
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
# fix up variable names in the data.frame.
# Variable names in R must start with a letter and contain letters or number or _.
# Blanks in variable names are not normally allowed. Blanks will be replaced by . (period)
cat("\nOriginal variable names in data frame\n")
##
## Original variable names in data frame
names(amph.df)
## [1] "file.name"
## [2] "Study Area Name"
## [3] "Transect Label"
## [4] "Date"
## [5] "Time"
## [6] "End Time"
## [7] "Surveyor"
## [8] "Species"
## [9] "Count"
## [10] "UTM Zone"
## [11] "Easting"
## [12] "Northing"
## [13] "Survey Type"
## [14] "Life Stage"
## [15] "Background Noise (for auditory surveys)"
## [16] "Calling Noise (for auditory surveys)"
## [17] "Preceeding 24 hour precipitation"
## [18] "Air temp (c)"
## [19] "Cloud cover"
## [20] "Wind speed"
## [21] "Snow coverage on ground"
## [22] "Current precipitation"
## [23] "BEC"
## [24] "Comments"
## [25] "Preceeding 24 hour precipitation (mm)"
names(amph.df) <- make.names(names(amph.df))
cat("\nCorrected variable names of data frame\n")
##
## Corrected variable names of data frame
names(amph.df)
## [1] "file.name"
## [2] "Study.Area.Name"
## [3] "Transect.Label"
## [4] "Date"
## [5] "Time"
## [6] "End.Time"
## [7] "Surveyor"
## [8] "Species"
## [9] "Count"
## [10] "UTM.Zone"
## [11] "Easting"
## [12] "Northing"
## [13] "Survey.Type"
## [14] "Life.Stage"
## [15] "Background.Noise..for.auditory.surveys."
## [16] "Calling.Noise..for.auditory.surveys."
## [17] "Preceeding.24.hour.precipitation"
## [18] "Air.temp..c."
## [19] "Cloud.cover"
## [20] "Wind.speed"
## [21] "Snow.coverage.on.ground"
## [22] "Current.precipitation"
## [23] "BEC"
## [24] "Comments"
## [25] "Preceeding.24.hour.precipitation..mm."
# Convert dates to R date format
xtabs(~Date, data=amph.df, exclude=NULL, na.action=na.pass) # check the date formats. Make sure that all yyyy-mm-dd
## Date
## 2013-05-22 2013-06-10 2014-04-09 2015-04-10
## 11 6 12 9
amph.df$Date <- as.Date(amph.df$Date, "%d-%b-%y", tz="UTC")
amph.df$Year <- as.numeric(format(amph.df$Date, "%Y"))
# Check that the Study Area Name is the same across all years
# Look at the output from the xtabs() to see if there are multiple spellings
# of the same Study.Area.Name.
# We will convert the Study.Area.Name to Proper Case.
amph.df$Study.Area.Name <- stringr::str_to_title(amph.df$Study.Area.Name)
xtabs(~Study.Area.Name, data=amph.df, exclude=NULL, na.action=na.pass)
## Study.Area.Name
## Alice Lake Park
## 38
xtabs(~Study.Area.Name+Year, data=amph.df, exclude=NULL, na.action=na.pass)
## Year
## Study.Area.Name 2013 2014 2015
## Alice Lake Park 17 12 9
# Check the Transect.Labels for typos
xtabs(~Study.Area.Name+Transect.Label, data=amph.df, exclude=NULL, na.action=na.pass)
## Transect.Label
## Study.Area.Name Fawn Lake Station Fawn Lake Transect
## Alice Lake Park 3 35
xtabs(~Transect.Label+Year, data=amph.df, exclude=NULL, na.action=na.pass)
## Year
## Transect.Label 2013 2014 2015
## Fawn Lake Station 2 1 0
## Fawn Lake Transect 15 11 9
# When is each station measured?
xtabs(~Transect.Label+Year, data=amph.df, exclude=NULL, na.action=na.pass)
## Year
## Transect.Label 2013 2014 2015
## Fawn Lake Station 2 1 0
## Fawn Lake Transect 15 11 9
xtabs(~Date+Transect.Label, data=amph.df, exclude=NULL, na.action=na.pass)
## Transect.Label
## Date Fawn Lake Station Fawn Lake Transect
## 2013-05-22 1 10
## 2013-06-10 1 5
## 2014-04-09 1 11
## 2015-04-10 0 9
# Check the Species code to make sure that they are all ok
# This isn't used anywhere in the analysis but is useful to know
xtabs(~Species, data=amph.df, exclude=NULL, na.action=na.pass)
## Species
## A-AMGR A-AMMA A-ANBO A-PRSE A-PSRE A-RAAU
## 5 2 6 8 9 5
## AMPHIBION Unidentifed
## 2 1
xtabs(~Species+Year, data=amph.df, exclude=NULL, na.action=na.pass)
## Year
## Species 2013 2014 2015
## A-AMGR 4 0 1
## A-AMMA 2 0 0
## A-ANBO 1 1 4
## A-PRSE 0 8 0
## A-PSRE 6 1 2
## A-RAAU 2 2 1
## AMPHIBION 2 0 0
## Unidentifed 0 0 1
# Check the survey type fiels
xtabs(~Year+Survey.Type, data=amph.df, exclude=NULL, na.action=na.pass)
## Survey.Type
## Year AU VI
## 2013 2 15
## 2014 1 11
## 2015 0 9
# Check the Life Stage Code
xtabs(~Year+Life.Stage, data=amph.df, exclude=NULL, na.action=na.pass)
## Life.Stage
## Year AD EG LA
## 2013 3 8 6
## 2014 3 9 0
## 2015 4 4 1
xtabs(~Survey.Type+Life.Stage, data=amph.df, exclude=NULL, na.action=na.pass)
## Life.Stage
## Survey.Type AD EG LA
## AU 3 0 0
## VI 7 21 7
# list the adult records
amph.df[ amph.df$Life.Stage=="AD",c("Study.Area.Name","Date","Survey.Type","Life.Stage","Count")]
## Study.Area.Name Date Survey.Type Life.Stage Count
## 1 Alice Lake Park 2013-05-22 AU AD 3
## 7 Alice Lake Park 2013-05-22 VI AD 2
## 17 Alice Lake Park 2013-06-10 AU AD 1
## 18 Alice Lake Park 2014-04-09 AU AD NA
## 19 Alice Lake Park 2014-04-09 VI AD 3
## 26 Alice Lake Park 2014-04-09 VI AD 1
## 33 Alice Lake Park 2015-04-10 VI AD 10
## 34 Alice Lake Park 2015-04-10 VI AD 5
## 36 Alice Lake Park 2015-04-10 VI AD 1
## 38 Alice Lake Park 2015-04-10 VI AD 1
# Get the file prefix
file.prefix <- make.names(amph.df$Study.Area.Name[1])
file.prefix <- gsub(".", '-', file.prefix, fixed=TRUE) # convert blanks to -
file.prefix <- file.path("Plots", file.prefix)
# Analysis of the calliing data
# Not possible at this time. See document.
# Analysis of the visual survey protocol
# Look at # egg masses
# We need to assume that detection probabilities are consistent over time and space
# There is only one monitoring station at Alice Lake so we don't need to model transect or year-specific
# process error. But I've generalized the code so that it also works with more than one station.
# Must impute 0 values for species not present in a year
dim(amph.df)
## [1] 38 26
amph.v.df <- amph.df[ amph.df$Survey.Type=='VI',]
dim(amph.v.df)
## [1] 35 26
# Count the total number of egg masses seen by species
eggmass.count.by.species <- plyr::ddply(amph.v.df, c("Study.Area.Name","Transect.Label","Year","Species"), summarize,
n.eggmass=sum(Life.Stage=='EG'))
eggmass.count.by.species
## Study.Area.Name Transect.Label Year Species n.eggmass
## 1 Alice Lake Park Fawn Lake Transect 2013 A-AMGR 3
## 2 Alice Lake Park Fawn Lake Transect 2013 A-AMMA 2
## 3 Alice Lake Park Fawn Lake Transect 2013 A-ANBO 0
## 4 Alice Lake Park Fawn Lake Transect 2013 A-PSRE 2
## 5 Alice Lake Park Fawn Lake Transect 2013 A-RAAU 1
## 6 Alice Lake Park Fawn Lake Transect 2013 AMPHIBION 0
## 7 Alice Lake Park Fawn Lake Transect 2014 A-ANBO 0
## 8 Alice Lake Park Fawn Lake Transect 2014 A-PRSE 7
## 9 Alice Lake Park Fawn Lake Transect 2014 A-RAAU 2
## 10 Alice Lake Park Fawn Lake Transect 2015 A-AMGR 1
## 11 Alice Lake Park Fawn Lake Transect 2015 A-ANBO 1
## 12 Alice Lake Park Fawn Lake Transect 2015 A-PSRE 1
## 13 Alice Lake Park Fawn Lake Transect 2015 A-RAAU 1
## 14 Alice Lake Park Fawn Lake Transect 2015 Unidentifed 0
# Count the total number of egg masses seen for ALL species
eggmass.count.all.species <- plyr::ddply(amph.v.df, c("Study.Area.Name","Transect.Label","Year","Species"), summarize,
n.eggmass=sum(Life.Stage=='EG'))
eggmass.count.all.species <- plyr::ddply(amph.v.df, c("Study.Area.Name","Transect.Label","Year"), summarize,
n.eggmass=sum(Life.Stage=='EG'))
eggmass.count.all.species$Species <- 'ALL.species'
eggmass.count <- rbind(eggmass.count.by.species, eggmass.count.all.species)
# Need to impute 0 values for species not seen in a year
# check out records where eggmass count not identified
select <- eggmass.count$Species %in% c("AMPHIBION","Unidentifed") # Notice type in unidentified
eggmass.count[ select,]
## Study.Area.Name Transect.Label Year Species n.eggmass
## 6 Alice Lake Park Fawn Lake Transect 2013 AMPHIBION 0
## 14 Alice Lake Park Fawn Lake Transect 2015 Unidentifed 0
eggmass.count <- eggmass.count[ !select, ]
# Create a record for each species x each year
complete.species.year <- expand.grid(Study.Area.Name=unique(eggmass.count$Study.Area.Name),
Transect.Label =unique(eggmass.count$Transect.Label),
Species=unique(eggmass.count$Species),
Year =unique(eggmass.count$Year, stringsAsFactors=FALSE))
eggmass.count <- merge(complete.species.year, eggmass.count, all.x=TRUE)
# replace all NA by 0
eggmass.count$n.eggmass[ is.na(eggmass.count$n.eggmass)] <- 0
# check that every species is given in each year
xtabs(~Species+Year, data=eggmass.count, exclude=NULL, na.action=na.pass)
## Year
## Species 2013 2014 2015
## A-AMGR 1 1 1
## A-AMMA 1 1 1
## A-ANBO 1 1 1
## A-PSRE 1 1 1
## A-RAAU 1 1 1
## A-PRSE 1 1 1
## ALL.species 1 1 1
# Section of code used for testing multiple transects#
#eggmass.count2 <- eggmass.count
#eggmass.count2$Transect.Label ='xx'
#eggmass.count2$n.eggmass <- rpois(1, eggmass.count$n.eggmass)
#eggmass.count2
#eggmass.count <- rbind(eggmass.count, eggmass.count2)
# Make a preliminary plot of total count by date
prelim.egg.plot <- ggplot(data=eggmass.count, aes(x=Year, y=n.eggmass, color=Transect.Label, linetype=Transect.Label))+
ggtitle("Amphibian egg mass count data")+
ylab("Amphibian egg masses")+
geom_point(position=position_dodge(width=.5))+
geom_line( position=position_dodge(width=.5))+
facet_wrap(~interaction(Study.Area.Name,Species), ncol=2, scales="free_y")+
scale_x_continuous(breaks=min(eggmass.count$Year,na.rm=TRUE):max(eggmass.count$Year,na.rm=TRUE))
prelim.egg.plot
ggplot2::ggsave(plot=prelim.egg.plot,
file=paste(file.prefix,'-plot-prelim-egg.png',sep=""),
h=6, w=6, units="in",dpi=300)
# This is a regression analysis with Year as the trend variable
# In this case, there is only one transect, so it is not necessary to have the Transect.Label in the model
# We also don't need to model process (year specific effects) because there is only 1 transect
# Fit a linear regression for each species
fits<- dlply(eggmass.count, "Species", function(eggmass.count){
cat("\n\n\n *** Starting analysis for species ", as.character(eggmass.count$Species[1]), "\n")
if(length(unique(eggmass.count$Transect.Label))>1){
eggmass.count$Transect.LabelF <- factor(eggmass.count$Transect.Label)
eggmass.count$YearF <- factor(eggmass.count$Year)
egg.fit <- lmerTest::lmer(log(n.eggmass) ~ Year + (1|Transect.LabelF) + (1|YearF) , data=eggmass.count)
print(anova(egg.fit, ddf='kenward-roger'))
print(summary(egg.fit))
print(VarCorr(egg.fit))
}
if(length(unique(eggmass.count$Transect.Label))==1){
# with only 1 transect, not necessary to put random effects effect in the model
egg.fit <- glm(n.eggmass ~ Year , data=eggmass.count, family=quasipoisson)
print(Anova(egg.fit, test="F", type=3))
print(summary(egg.fit))
}
list(Species=eggmass.count$Species[1], fit=egg.fit)
})
##
##
##
## *** Starting analysis for species A-AMGR
## Analysis of Deviance Table (Type III tests)
##
## Response: n.eggmass
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 1.5790 1 0.8759 0.5211
## Residuals 1.8028 1
##
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
##
## Deviance Residuals:
## 1 2 3
## 0.3296 -1.4631 0.6796
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1679.9760 1946.1905 0.863 0.547
## Year -0.8341 0.9666 -0.863 0.547
##
## (Dispersion parameter for quasipoisson family taken to be 1.802804)
##
## Null deviance: 4.2902 on 2 degrees of freedom
## Residual deviance: 2.7112 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
##
##
##
##
## *** Starting analysis for species A-AMMA
## Analysis of Deviance Table (Type III tests)
##
## Response: n.eggmass
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 4.3944 1 2.5127e+10 4.016e-06 ***
## Residuals 0.0000 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
##
## Deviance Residuals:
## 1 2 3
## 2.11e-08 -1.87e-05 -2.11e-08
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 46621.82 2013.00 23.16 0.0275 *
## Year -23.16 1.00 -23.16 0.0275 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 4.753932e-10)
##
## Null deviance: 4.3944e+00 on 2 degrees of freedom
## Residual deviance: 3.4978e-10 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 22
##
##
##
##
## *** Starting analysis for species A-ANBO
## Analysis of Deviance Table (Type III tests)
##
## Response: n.eggmass
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 2.1972 1 1.4968e+10 5.204e-06 ***
## Residuals 0.0000 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
##
## Deviance Residuals:
## 1 2 3
## -2.110e-08 -1.713e-05 0.000e+00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -45623.57 2015.00 -22.64 0.0281 *
## Year 22.64 1.00 22.64 0.0281 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 3.990352e-10)
##
## Null deviance: 2.1972e+00 on 2 degrees of freedom
## Residual deviance: 2.9359e-10 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 22
##
##
##
##
## *** Starting analysis for species A-PSRE
## Analysis of Deviance Table (Type III tests)
##
## Response: n.eggmass
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 0.51094 1 0.3558 0.6576
## Residuals 1.43614 1
##
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
##
## Deviance Residuals:
## 1 2 3
## 0.352 -1.353 0.555
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1052.1098 1823.2336 0.577 0.667
## Year -0.5224 0.9054 -0.577 0.667
##
## (Dispersion parameter for quasipoisson family taken to be 1.436144)
##
## Null deviance: 2.7726 on 2 degrees of freedom
## Residual deviance: 2.2616 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
##
##
##
##
## *** Starting analysis for species A-RAAU
## Analysis of Deviance Table (Type III tests)
##
## Response: n.eggmass
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 0.0 1 0 1
## Residuals 0.5 1
##
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
##
## Deviance Residuals:
## 1 2 3
## -0.3022 0.5372 -0.3022
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.877e-01 8.721e+02 0 1
## Year 1.050e-13 4.330e-01 0 1
##
## (Dispersion parameter for quasipoisson family taken to be 0.5000005)
##
## Null deviance: 0.47113 on 2 degrees of freedom
## Residual deviance: 0.47113 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
##
##
##
##
## *** Starting analysis for species A-PRSE
## Analysis of Deviance Table (Type III tests)
##
## Response: n.eggmass
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 0 1 0 1
## Residuals 14 1
##
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
##
## Deviance Residuals:
## 1 2 3
## -2.160 2.459 -2.160
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.473e-01 3.488e+03 0 1
## Year -4.488e-13 1.732e+00 0 1
##
## (Dispersion parameter for quasipoisson family taken to be 14.00001)
##
## Null deviance: 15.381 on 2 degrees of freedom
## Residual deviance: 15.381 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 6
##
##
##
##
## *** Starting analysis for species ALL.species
## Analysis of Deviance Table (Type III tests)
##
## Response: n.eggmass
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 1.1507 1 1.0721 0.4889
## Residuals 1.0733 1
##
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
##
## Deviance Residuals:
## 1 2 3
## -0.3711 0.8001 -0.5046
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 585.3617 569.3084 1.028 0.491
## Year -0.2897 0.2827 -1.025 0.492
##
## (Dispersion parameter for quasipoisson family taken to be 1.073323)
##
## Null deviance: 2.1832 on 2 degrees of freedom
## Residual deviance: 1.0325 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
# Get the model diagnostic plots
l_ply(fits, function(x){
cat("\n\n\n *** Diagnostic plot for species ", as.character(x$Species), "\n")
if(length(unique(x$fit$data$Transect.Label))>1){
diag.plot <- sf.autoplot.lmer(x$fit) # residual and other diagnostic plots
plot(diag.plot)
}
if(length(unique(x$fit$data$Transect.Label))==1){
diag.plot <- autoplot(x$fit) # residual and other diagnostic plots
show(diag.plot)
}
ggplot2:: ggsave(plot=diag.plot,
file=paste(file.prefix,"-egg-residual-plot-",x$Species,".png",sep=""),
h=6, w=6, units="in", dpi=300)
})
##
##
##
## *** Diagnostic plot for species A-AMGR
##
##
##
## *** Diagnostic plot for species A-AMMA
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
##
##
##
## *** Diagnostic plot for species A-ANBO
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
##
##
##
## *** Diagnostic plot for species A-PSRE
##
##
##
## *** Diagnostic plot for species A-RAAU
##
##
##
## *** Diagnostic plot for species A-PRSE
##
##
##
## *** Diagnostic plot for species ALL.species
l_ply(fits, function(x){
cat("\n\n\n *** Check for autocorrelation for species ", as.character(x$Species), "\n")
eggmass.count <- x$fit$data
# check for autocorrelation
if(length(unique(eggmass.count$Transect.Label))>1){
eggmass.count$resid <- log(eggmass.count$n.eggmass) - predict(x$fit, newdata=eggmass.count, re.form=~0)
}
#browser()
if(length(unique(eggmass.count$Transect.Label))==1){
eggmass.count$resid <- log(eggmass.count$n.eggmass+.1) - predict(x$fit, newdata=eggmass.count, type="link")
}
mean.resid <- plyr::ddply(eggmass.count, "Year", summarize, mean.resid=mean(resid))
resid.fit <- lm( mean.resid ~ 1, data=mean.resid)
dwres1 <- car::durbinWatsonTest(resid.fit)
print(dwres1)
dwres2 <- lmtest::dwtest(resid.fit)
print(dwres2)
})
##
##
##
## *** Check for autocorrelation for species A-AMGR
## lag Autocorrelation D-W Statistic p-value
## 1 -0.6439684 2.931905 0.518
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 2.9319, p-value = 0.8818
## alternative hypothesis: true autocorrelation is greater than 0
##
##
##
##
## *** Check for autocorrelation for species A-AMMA
## lag Autocorrelation D-W Statistic p-value
## 1 -0.001098058 1.003294 0.122
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 1.0033, p-value = 0.01854
## alternative hypothesis: true autocorrelation is greater than 0
##
##
##
##
## *** Check for autocorrelation for species A-ANBO
## lag Autocorrelation D-W Statistic p-value
## 1 -0.0006940065 1.002082 NA
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 1.0021, p-value = 0.01282
## alternative hypothesis: true autocorrelation is greater than 0
##
##
##
##
## *** Check for autocorrelation for species A-PSRE
## lag Autocorrelation D-W Statistic p-value
## 1 -0.6561264 2.968379 0.432
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 2.9684, p-value = 0.9197
## alternative hypothesis: true autocorrelation is greater than 0
##
##
##
##
## *** Check for autocorrelation for species A-RAAU
## lag Autocorrelation D-W Statistic p-value
## 1 -0.6666667 3 0.166
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 3, p-value = 1
## alternative hypothesis: true autocorrelation is greater than 0
##
##
##
##
## *** Check for autocorrelation for species A-PRSE
## lag Autocorrelation D-W Statistic p-value
## 1 -0.6666667 3 0
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 3, p-value = 1
## alternative hypothesis: true autocorrelation is greater than 0
##
##
##
##
## *** Check for autocorrelation for species ALL.species
## lag Autocorrelation D-W Statistic p-value
## 1 -0.642872 2.928616 0.444
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 2.9286, p-value = 0.8789
## alternative hypothesis: true autocorrelation is greater than 0
# extract a table of statistics for each study area x species combiations
egg.slopes <- ldply(fits, function(x){
eggmass.count <- x$fit$data
egg.fit <- x$fit
if(length(unique(eggmass.count$Transect.Label))==1){
egg.slopes <- data.frame(
Study.Area.Name =eggmass.count$Study.Area.Name[1],
slope = coef(egg.fit)["Year"],
slope.se = sqrt(diag(vcov(egg.fit)))["Year"],
p.value = summary(egg.fit)$coefficients[row.names(summary(egg.fit)$coefficients)=="Year" ,"Pr(>|t|)"],
stringsAsFactors=FALSE
)
}
if(length(unique(eggmass.count$Transect.Label))>1){
egg.slopes <- data.frame(
Study.Area.Name =eggmass.count$Study.Area.Name[1],
slope = fixef(egg.fit)["Year"],
slope.se = summary(egg.fit)$coefficients["Year","Pr(>|t|)"],
p.value = summary(egg.fit)$coefficients[row.names(summary(egg.fit)$coefficients)=="Year" ,"Pr(>|t|)"],
stringsAsFactors=FALSE)
}
egg.slopes
})
egg.slopes
## Species Study.Area.Name slope slope.se p.value
## 1 A-AMGR Alice Lake Park -8.341152e-01 0.9665708 0.54674457
## 2 A-AMMA Alice Lake Park -2.316003e+01 0.9999993 0.02747079
## 3 A-ANBO Alice Lake Park 2.264197e+01 0.9999992 0.02809852
## 4 A-PSRE Alice Lake Park -5.224423e-01 0.9054296 0.66682897
## 5 A-RAAU Alice Lake Park 1.049719e-13 0.4330127 1.00000000
## 6 A-PRSE Alice Lake Park -4.487664e-13 1.7320508 1.00000000
## 7 ALL.species Alice Lake Park -2.896940e-01 0.2827022 0.49222412
# compute the fitted values from the model
egg.fitted <- ldply(fits, function(x){
eggmass.count <- x$fit$data
egg.fit <- x$fit
Transect.Label <- unique(as.character(eggmass.count$Transect.Label))
newdata <- expand.grid(Year=seq(min(eggmass.count$Year, na.rm=TRUE),max(eggmass.count$Year, na.rm=TRUE), .1),
Transect.LabelF=Transect.Label)
newdata$Transect.Label <-as.character(newdata$Transect.LabelF)
if(length(unique(eggmass.count$Transect.Label))==1){
newdata$pred.mean <- predict(egg.fit, newdata=newdata,type="response")
}
if(length(unique(eggmass.count$Transect.Label))>1){
newdata$pred.mean <- exp(predict(egg.fit, newdata=newdata,type="response", re.form=~0))
}
# we now must average over all of the transect labels
egg.fitted <- plyr::ddply(newdata, c("Year","Transect.Label"), plyr::summarize, pred.mean=mean(pred.mean))
egg.fitted$Study.Area.Name <- eggmass.count$Study.Area.Name
egg.fitted
})
egg.fitted
## Species Year Transect.Label pred.mean Study.Area.Name
## 1 A-AMGR 2013.0 Fawn Lake Transect 2.464816e+00 Alice Lake Park
## 2 A-AMGR 2013.1 Fawn Lake Transect 2.267563e+00 Alice Lake Park
## 3 A-AMGR 2013.2 Fawn Lake Transect 2.086096e+00 Alice Lake Park
## 4 A-AMGR 2013.3 Fawn Lake Transect 1.919151e+00 Alice Lake Park
## 5 A-AMGR 2013.4 Fawn Lake Transect 1.765566e+00 Alice Lake Park
## 6 A-AMGR 2013.5 Fawn Lake Transect 1.624272e+00 Alice Lake Park
## 7 A-AMGR 2013.6 Fawn Lake Transect 1.494285e+00 Alice Lake Park
## 8 A-AMGR 2013.7 Fawn Lake Transect 1.374702e+00 Alice Lake Park
## 9 A-AMGR 2013.8 Fawn Lake Transect 1.264688e+00 Alice Lake Park
## 10 A-AMGR 2013.9 Fawn Lake Transect 1.163478e+00 Alice Lake Park
## 11 A-AMGR 2014.0 Fawn Lake Transect 1.070368e+00 Alice Lake Park
## 12 A-AMGR 2014.1 Fawn Lake Transect 9.847087e-01 Alice Lake Park
## 13 A-AMGR 2014.2 Fawn Lake Transect 9.059049e-01 Alice Lake Park
## 14 A-AMGR 2014.3 Fawn Lake Transect 8.334076e-01 Alice Lake Park
## 15 A-AMGR 2014.4 Fawn Lake Transect 7.667120e-01 Alice Lake Park
## 16 A-AMGR 2014.5 Fawn Lake Transect 7.053540e-01 Alice Lake Park
## 17 A-AMGR 2014.6 Fawn Lake Transect 6.489062e-01 Alice Lake Park
## 18 A-AMGR 2014.7 Fawn Lake Transect 5.969759e-01 Alice Lake Park
## 19 A-AMGR 2014.8 Fawn Lake Transect 5.492014e-01 Alice Lake Park
## 20 A-AMGR 2014.9 Fawn Lake Transect 5.052502e-01 Alice Lake Park
## 21 A-AMGR 2015.0 Fawn Lake Transect 4.648162e-01 Alice Lake Park
## 22 A-AMMA 2013.0 Fawn Lake Transect 2.000000e+00 Alice Lake Park
## 23 A-AMMA 2013.1 Fawn Lake Transect 1.973344e-01 Alice Lake Park
## 24 A-AMMA 2013.2 Fawn Lake Transect 1.947044e-02 Alice Lake Park
## 25 A-AMMA 2013.3 Fawn Lake Transect 1.921094e-03 Alice Lake Park
## 26 A-AMMA 2013.4 Fawn Lake Transect 1.895490e-04 Alice Lake Park
## 27 A-AMMA 2013.5 Fawn Lake Transect 1.870227e-05 Alice Lake Park
## 28 A-AMMA 2013.6 Fawn Lake Transect 1.845301e-06 Alice Lake Park
## 29 A-AMMA 2013.7 Fawn Lake Transect 1.820707e-07 Alice Lake Park
## 30 A-AMMA 2013.8 Fawn Lake Transect 1.796441e-08 Alice Lake Park
## 31 A-AMMA 2013.9 Fawn Lake Transect 1.772498e-09 Alice Lake Park
## 32 A-AMMA 2014.0 Fawn Lake Transect 1.748875e-10 Alice Lake Park
## 33 A-AMMA 2014.1 Fawn Lake Transect 1.725566e-11 Alice Lake Park
## 34 A-AMMA 2014.2 Fawn Lake Transect 1.702568e-12 Alice Lake Park
## 35 A-AMMA 2014.3 Fawn Lake Transect 1.679876e-13 Alice Lake Park
## 36 A-AMMA 2014.4 Fawn Lake Transect 1.657487e-14 Alice Lake Park
## 37 A-AMMA 2014.5 Fawn Lake Transect 1.635397e-15 Alice Lake Park
## 38 A-AMMA 2014.6 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 39 A-AMMA 2014.7 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 40 A-AMMA 2014.8 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 41 A-AMMA 2014.9 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 42 A-AMMA 2015.0 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 43 A-ANBO 2013.0 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 44 A-ANBO 2013.1 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 45 A-ANBO 2013.2 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 46 A-ANBO 2013.3 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 47 A-ANBO 2013.4 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 48 A-ANBO 2013.5 Fawn Lake Transect 1.778588e-15 Alice Lake Park
## 49 A-ANBO 2013.6 Fawn Lake Transect 1.711606e-14 Alice Lake Park
## 50 A-ANBO 2013.7 Fawn Lake Transect 1.647146e-13 Alice Lake Park
## 51 A-ANBO 2013.8 Fawn Lake Transect 1.585113e-12 Alice Lake Park
## 52 A-ANBO 2013.9 Fawn Lake Transect 1.525417e-11 Alice Lake Park
## 53 A-ANBO 2014.0 Fawn Lake Transect 1.467969e-10 Alice Lake Park
## 54 A-ANBO 2014.1 Fawn Lake Transect 1.412685e-09 Alice Lake Park
## 55 A-ANBO 2014.2 Fawn Lake Transect 1.359482e-08 Alice Lake Park
## 56 A-ANBO 2014.3 Fawn Lake Transect 1.308283e-07 Alice Lake Park
## 57 A-ANBO 2014.4 Fawn Lake Transect 1.259013e-06 Alice Lake Park
## 58 A-ANBO 2014.5 Fawn Lake Transect 1.211598e-05 Alice Lake Park
## 59 A-ANBO 2014.6 Fawn Lake Transect 1.165968e-04 Alice Lake Park
## 60 A-ANBO 2014.7 Fawn Lake Transect 1.122057e-03 Alice Lake Park
## 61 A-ANBO 2014.8 Fawn Lake Transect 1.079800e-02 Alice Lake Park
## 62 A-ANBO 2014.9 Fawn Lake Transect 1.039134e-01 Alice Lake Park
## 63 A-ANBO 2015.0 Fawn Lake Transect 1.000000e+00 Alice Lake Park
## 64 A-PSRE 2013.0 Fawn Lake Transect 1.542573e+00 Alice Lake Park
## 65 A-PSRE 2013.1 Fawn Lake Transect 1.464051e+00 Alice Lake Park
## 66 A-PSRE 2013.2 Fawn Lake Transect 1.389527e+00 Alice Lake Park
## 67 A-PSRE 2013.3 Fawn Lake Transect 1.318796e+00 Alice Lake Park
## 68 A-PSRE 2013.4 Fawn Lake Transect 1.251665e+00 Alice Lake Park
## 69 A-PSRE 2013.5 Fawn Lake Transect 1.187952e+00 Alice Lake Park
## 70 A-PSRE 2013.6 Fawn Lake Transect 1.127481e+00 Alice Lake Park
## 71 A-PSRE 2013.7 Fawn Lake Transect 1.070089e+00 Alice Lake Park
## 72 A-PSRE 2013.8 Fawn Lake Transect 1.015619e+00 Alice Lake Park
## 73 A-PSRE 2013.9 Fawn Lake Transect 9.639206e-01 Alice Lake Park
## 74 A-PSRE 2014.0 Fawn Lake Transect 9.148542e-01 Alice Lake Park
## 75 A-PSRE 2014.1 Fawn Lake Transect 8.682854e-01 Alice Lake Park
## 76 A-PSRE 2014.2 Fawn Lake Transect 8.240871e-01 Alice Lake Park
## 77 A-PSRE 2014.3 Fawn Lake Transect 7.821387e-01 Alice Lake Park
## 78 A-PSRE 2014.4 Fawn Lake Transect 7.423255e-01 Alice Lake Park
## 79 A-PSRE 2014.5 Fawn Lake Transect 7.045389e-01 Alice Lake Park
## 80 A-PSRE 2014.6 Fawn Lake Transect 6.686758e-01 Alice Lake Park
## 81 A-PSRE 2014.7 Fawn Lake Transect 6.346382e-01 Alice Lake Park
## 82 A-PSRE 2014.8 Fawn Lake Transect 6.023333e-01 Alice Lake Park
## 83 A-PSRE 2014.9 Fawn Lake Transect 5.716727e-01 Alice Lake Park
## 84 A-PSRE 2015.0 Fawn Lake Transect 5.425729e-01 Alice Lake Park
## 85 A-RAAU 2013.0 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 86 A-RAAU 2013.1 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 87 A-RAAU 2013.2 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 88 A-RAAU 2013.3 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 89 A-RAAU 2013.4 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 90 A-RAAU 2013.5 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 91 A-RAAU 2013.6 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 92 A-RAAU 2013.7 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 93 A-RAAU 2013.8 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 94 A-RAAU 2013.9 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 95 A-RAAU 2014.0 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 96 A-RAAU 2014.1 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 97 A-RAAU 2014.2 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 98 A-RAAU 2014.3 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 99 A-RAAU 2014.4 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 100 A-RAAU 2014.5 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 101 A-RAAU 2014.6 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 102 A-RAAU 2014.7 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 103 A-RAAU 2014.8 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 104 A-RAAU 2014.9 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 105 A-RAAU 2015.0 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 106 A-PRSE 2013.0 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 107 A-PRSE 2013.1 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 108 A-PRSE 2013.2 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 109 A-PRSE 2013.3 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 110 A-PRSE 2013.4 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 111 A-PRSE 2013.5 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 112 A-PRSE 2013.6 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 113 A-PRSE 2013.7 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 114 A-PRSE 2013.8 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 115 A-PRSE 2013.9 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 116 A-PRSE 2014.0 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 117 A-PRSE 2014.1 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 118 A-PRSE 2014.2 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 119 A-PRSE 2014.3 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 120 A-PRSE 2014.4 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 121 A-PRSE 2014.5 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 122 A-PRSE 2014.6 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 123 A-PRSE 2014.7 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 124 A-PRSE 2014.8 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 125 A-PRSE 2014.9 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 126 A-PRSE 2015.0 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 127 ALL.species 2013.0 Fawn Lake Transect 9.095895e+00 Alice Lake Park
## 128 ALL.species 2013.1 Fawn Lake Transect 8.836173e+00 Alice Lake Park
## 129 ALL.species 2013.2 Fawn Lake Transect 8.583866e+00 Alice Lake Park
## 130 ALL.species 2013.3 Fawn Lake Transect 8.338764e+00 Alice Lake Park
## 131 ALL.species 2013.4 Fawn Lake Transect 8.100661e+00 Alice Lake Park
## 132 ALL.species 2013.5 Fawn Lake Transect 7.869356e+00 Alice Lake Park
## 133 ALL.species 2013.6 Fawn Lake Transect 7.644656e+00 Alice Lake Park
## 134 ALL.species 2013.7 Fawn Lake Transect 7.426372e+00 Alice Lake Park
## 135 ALL.species 2013.8 Fawn Lake Transect 7.214321e+00 Alice Lake Park
## 136 ALL.species 2013.9 Fawn Lake Transect 7.008324e+00 Alice Lake Park
## 137 ALL.species 2014.0 Fawn Lake Transect 6.808210e+00 Alice Lake Park
## 138 ALL.species 2014.1 Fawn Lake Transect 6.613810e+00 Alice Lake Park
## 139 ALL.species 2014.2 Fawn Lake Transect 6.424960e+00 Alice Lake Park
## 140 ALL.species 2014.3 Fawn Lake Transect 6.241503e+00 Alice Lake Park
## 141 ALL.species 2014.4 Fawn Lake Transect 6.063285e+00 Alice Lake Park
## 142 ALL.species 2014.5 Fawn Lake Transect 5.890155e+00 Alice Lake Park
## 143 ALL.species 2014.6 Fawn Lake Transect 5.721968e+00 Alice Lake Park
## 144 ALL.species 2014.7 Fawn Lake Transect 5.558584e+00 Alice Lake Park
## 145 ALL.species 2014.8 Fawn Lake Transect 5.399866e+00 Alice Lake Park
## 146 ALL.species 2014.9 Fawn Lake Transect 5.245679e+00 Alice Lake Park
## 147 ALL.species 2015.0 Fawn Lake Transect 5.095895e+00 Alice Lake Park
# Plot with trend line
egg.plot.summary <- ggplot2::ggplot(data=eggmass.count,
aes(x=Year, y=n.eggmass))+
ggtitle("Amphibian eggmass count ")+
ylab("Amphibian Eggmass Count")+
geom_point(size=3, aes(color=Transect.Label))+
facet_wrap(~interaction(Study.Area.Name,Species), ncol=2, scales="free_y")+
geom_line(data=egg.fitted, aes(x=Year,y=pred.mean))+
scale_x_continuous(breaks=min(eggmass.count$Year,na.rm=TRUE):max(eggmass.count$Year,na.rm=TRUE))+
geom_text(data=egg.slopes, aes(x=min(eggmass.count$Year, na.rm=TRUE), y=max(eggmass.count$n.eggmass, na.rm=TRUE)),
label=paste("Slope (on log scale) : ",round(egg.slopes$slope,2),
" ( SE " ,round(egg.slopes$slope.se,2),")",
" \np :" ,round(egg.slopes$p.value,3)),
hjust="left",vjust="top")+
theme(legend.position="top")
egg.plot.summary
ggsave(plot=egg.plot.summary,
file=paste(file.prefix,'-egg-plot-summary.png',sep=""),
h=6, w=6, units="in", dpi=300)
# Analysis of the visual survey protocol
# Look at # adults
# We need to assume that detection probabilities are consistent over time and space
# There is only one monitoring station at Alice Lake
dim(amph.df)
## [1] 38 26
amph.v.df <- amph.df[ amph.df$Survey.Type=='VI',]
dim(amph.v.df)
## [1] 35 26
# Count the total number of egg masses seen by species
adults.count.by.species <- plyr::ddply(amph.v.df, c("Study.Area.Name","Transect.Label","Year","Species"), summarize,
n.adults=sum(Life.Stage=='AD'))
adults.count.by.species
## Study.Area.Name Transect.Label Year Species n.adults
## 1 Alice Lake Park Fawn Lake Transect 2013 A-AMGR 0
## 2 Alice Lake Park Fawn Lake Transect 2013 A-AMMA 0
## 3 Alice Lake Park Fawn Lake Transect 2013 A-ANBO 0
## 4 Alice Lake Park Fawn Lake Transect 2013 A-PSRE 0
## 5 Alice Lake Park Fawn Lake Transect 2013 A-RAAU 1
## 6 Alice Lake Park Fawn Lake Transect 2013 AMPHIBION 0
## 7 Alice Lake Park Fawn Lake Transect 2014 A-ANBO 1
## 8 Alice Lake Park Fawn Lake Transect 2014 A-PRSE 1
## 9 Alice Lake Park Fawn Lake Transect 2014 A-RAAU 0
## 10 Alice Lake Park Fawn Lake Transect 2015 A-AMGR 0
## 11 Alice Lake Park Fawn Lake Transect 2015 A-ANBO 3
## 12 Alice Lake Park Fawn Lake Transect 2015 A-PSRE 1
## 13 Alice Lake Park Fawn Lake Transect 2015 A-RAAU 0
## 14 Alice Lake Park Fawn Lake Transect 2015 Unidentifed 0
# Count the total number of egg masses seen for ALL species
adults.count.all.species <- plyr::ddply(amph.v.df, c("Study.Area.Name","Transect.Label","Year","Species"), summarize,
n.adults=sum(Life.Stage=='AD'))
adults.count.all.species <- plyr::ddply(amph.v.df, c("Study.Area.Name","Transect.Label","Year"), summarize,
n.adults=sum(Life.Stage=='AD'))
adults.count.all.species$Species <- 'ALL.species'
adults.count <- rbind(adults.count.by.species, adults.count.all.species)
# Need to impute 0 values for species not seen in a year
# check out records where adults count not identified
select <- adults.count$Species %in% c("AMPHIBION","Unidentifed") # Notice type in unidentified
adults.count[ select,]
## Study.Area.Name Transect.Label Year Species n.adults
## 6 Alice Lake Park Fawn Lake Transect 2013 AMPHIBION 0
## 14 Alice Lake Park Fawn Lake Transect 2015 Unidentifed 0
adults.count <- adults.count[ !select, ]
# Create a record for each species x each year
complete.species.year <- expand.grid(Study.Area.Name=unique(adults.count$Study.Area.Name),
Transect.Label =unique(adults.count$Transect.Label),
Species=unique(adults.count$Species),
Year =unique(adults.count$Year, stringsAsFactors=FALSE))
adults.count <- merge(complete.species.year, adults.count, all.x=TRUE)
# replace all NA by 0
adults.count$n.adults[ is.na(adults.count$n.adults)] <- 0
# check that every species is given in each year
xtabs(~Species+Year, data=adults.count, exclude=NULL, na.action=na.pass)
## Year
## Species 2013 2014 2015
## A-AMGR 1 1 1
## A-AMMA 1 1 1
## A-ANBO 1 1 1
## A-PSRE 1 1 1
## A-RAAU 1 1 1
## A-PRSE 1 1 1
## ALL.species 1 1 1
# Section of code used for testing multiple transects#
#adults.count2 <- adults.count
#adults.count2$Transect.Label ='xx'
#adults.count2$n.adults <- rpois(1, 1+adults.count$n.adults)
#adults.count2
#adults.count <- rbind(adults.count, adults.count2)
# Make a preliminary plot of total count by date
# Make a preliminary plot of total count by date
prelim.adult.plot <- ggplot(data=adults.count, aes(x=Year, y=n.adults, color=Transect.Label, linetype=Transect.Label))+
ggtitle("Amphibian adult count data")+
ylab("Amphibian adults")+
geom_point(position=position_dodge(width=.5))+
geom_line( position=position_dodge(width=.5))+
facet_wrap(~interaction(Study.Area.Name,Species), ncol=2, scales="free_y")+
scale_x_continuous(breaks=min(adults.count$Year,na.rm=TRUE):max(adults.count$Year,na.rm=TRUE))
prelim.adult.plot
ggplot2::ggsave(plot=prelim.adult.plot,
file=paste(file.prefix,'-plot-prelim-adult.png',sep=""),
h=6, w=6, units="in",dpi=300)
# This is a regression analysis with Year as the trend variable
# In this case, there is only one transect, so it is not necessary to have the Transect.Label in the model
# We also don't need to model process (year specific effects) because there is only 1 transect
# Fit a linear regression for each species
fits<- dlply(adults.count, "Species", function(adults.count){
cat("\n\n\n *** Starting analysis for species ", as.character(adults.count$Species[1]), "\n")
if(length(unique(adults.count$Transect.Label))>1){
adults.count$Transect.LabelF <- factor(adults.count$Transect.Label)
adults.count$YearF <- factor(adults.count$Year)
adult.fit <- lmerTest::lmer(log(n.adults) ~ Year + (1|Transect.LabelF) + (1|YearF) , data=adults.count)
print(anova(adult.fit, ddf='kenward-roger'))
print(summary(adult.fit))
print(VarCorr(adult.fit))
}
if(length(unique(adults.count$Transect.Label))==1){
# with only 1 transect, not necessary to put random effects effect in the model
adult.fit <- glm(n.adults ~ Year , data=adults.count, family=quasipoisson)
print(Anova(adult.fit, test="F", type=3))
print(summary(adult.fit))
}
list(Species=adults.count$Species[1], fit=adult.fit)
})
##
##
##
## *** Starting analysis for species A-AMGR
## Analysis of Deviance Table (Type III tests)
##
## Response: n.adults
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 0.0000e+00 1 0 1
## Residuals 2.2748e-10 1
##
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
##
## Deviance Residuals:
## 1 2 3
## -1.231e-05 -1.231e-05 -1.231e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -23.303 2466.636 -0.009 0.994
## Year 0.000 1.225 0.000 1.000
##
## (Dispersion parameter for quasipoisson family taken to be 6.183461e-10)
##
## Null deviance: 0.0000e+00 on 2 degrees of freedom
## Residual deviance: 4.5495e-10 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 21
##
##
##
##
## *** Starting analysis for species A-AMMA
## Analysis of Deviance Table (Type III tests)
##
## Response: n.adults
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 0.0000e+00 1 0 1
## Residuals 2.2748e-10 1
##
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
##
## Deviance Residuals:
## 1 2 3
## -1.231e-05 -1.231e-05 -1.231e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -23.303 2466.636 -0.009 0.994
## Year 0.000 1.225 0.000 1.000
##
## (Dispersion parameter for quasipoisson family taken to be 6.183461e-10)
##
## Null deviance: 0.0000e+00 on 2 degrees of freedom
## Residual deviance: 4.5495e-10 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 21
##
##
##
##
## *** Starting analysis for species A-ANBO
## Analysis of Deviance Table (Type III tests)
##
## Response: n.adults
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 3.8586 1 13.083 0.1717
## Residuals 0.2949 1
##
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
##
## Deviance Residuals:
## 1 2 3
## -0.55294 0.34401 -0.08681
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3048.0126 1064.8901 -2.862 0.214
## Year 1.5132 0.5285 2.863 0.214
##
## (Dispersion parameter for quasipoisson family taken to be 0.294941)
##
## Null deviance: 4.29022 on 2 degrees of freedom
## Residual deviance: 0.43162 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
##
##
##
##
## *** Starting analysis for species A-PSRE
## Analysis of Deviance Table (Type III tests)
##
## Response: n.adults
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 2.1972 1 1.4968e+10 5.204e-06 ***
## Residuals 0.0000 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
##
## Deviance Residuals:
## 1 2 3
## -2.110e-08 -1.713e-05 0.000e+00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -45623.57 2015.00 -22.64 0.0281 *
## Year 22.64 1.00 22.64 0.0281 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 3.990352e-10)
##
## Null deviance: 2.1972e+00 on 2 degrees of freedom
## Residual deviance: 2.9359e-10 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 22
##
##
##
##
## *** Starting analysis for species A-RAAU
## Analysis of Deviance Table (Type III tests)
##
## Response: n.adults
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 2.1972 1 1.4968e+10 5.204e-06 ***
## Residuals 0.0000 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
##
## Deviance Residuals:
## 1 2 3
## 0.000e+00 -1.713e-05 -2.110e-08
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45578.29 2013.00 22.64 0.0281 *
## Year -22.64 1.00 -22.64 0.0281 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 3.990352e-10)
##
## Null deviance: 2.1972e+00 on 2 degrees of freedom
## Residual deviance: 2.9359e-10 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 22
##
##
##
##
## *** Starting analysis for species A-PRSE
## Analysis of Deviance Table (Type III tests)
##
## Response: n.adults
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 0 1 0 1
## Residuals 2 1
##
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
##
## Deviance Residuals:
## 1 2 3
## -0.8165 0.9295 -0.8165
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.099e+00 3.488e+03 0 1
## Year 4.340e-13 1.732e+00 0 1
##
## (Dispersion parameter for quasipoisson family taken to be 2.000103)
##
## Null deviance: 2.1972 on 2 degrees of freedom
## Residual deviance: 2.1972 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
##
##
##
##
## *** Starting analysis for species ALL.species
## Analysis of Deviance Table (Type III tests)
##
## Response: n.adults
## Error estimate based on Pearson residuals
##
## SS Df F Pr(>F)
## Year 2.0008 1 6.6857e+21 7.786e-12 ***
## Residuals 0.0000 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
##
## Deviance Residuals:
## 1 2 3
## 0 0 0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.395e+03 1.808e-08 -7.717e+10 8.25e-12 ***
## Year 6.931e-01 8.976e-12 7.722e+10 8.24e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 2.992646e-22)
##
## Null deviance: 2.0008e+00 on 2 degrees of freedom
## Residual deviance: -2.9926e-22 on 1 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 3
# Get the model diagnostic plots
l_ply(fits, function(x){
cat("\n\n\n *** Diagnostic plot for species ", as.character(x$Species), "\n")
if(length(unique(x$fit$data$Transect.Label))>1){
diag.plot <- sf.autoplot.lmer(x$fit) # residual and other diagnostic plots
plot(diag.plot)
}
if(length(unique(x$fit$data$Transect.Label))==1){
diag.plot <- autoplot(x$fit) # residual and other diagnostic plots
show(diag.plot)
}
ggplot2:: ggsave(plot=diag.plot,
file=paste(file.prefix,"-adult-residual-plot-",x$Species,".png",sep=""),
h=6, w=6, units="in", dpi=300)
})
##
##
##
## *** Diagnostic plot for species A-AMGR
##
##
##
## *** Diagnostic plot for species A-AMMA
##
##
##
## *** Diagnostic plot for species A-ANBO
##
##
##
## *** Diagnostic plot for species A-PSRE
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
##
##
##
## *** Diagnostic plot for species A-RAAU
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
##
##
##
## *** Diagnostic plot for species A-PRSE
##
##
##
## *** Diagnostic plot for species ALL.species
l_ply(fits, function(x){
cat("\n\n\n *** Check for autocorrelation for species ", as.character(x$Species), "\n")
adults.count <- x$fit$data
# check for autocorrelation
if(length(unique(adults.count$Transect.Label))>1){
adults.count$resid <- log(adults.count$n.adults) - predict(x$fit, newdata=adults.count, re.form=~0)
}
#browser()
if(length(unique(adults.count$Transect.Label))==1){
adults.count$resid <- log(adults.count$n.adults+.1) - predict(x$fit, newdata=adults.count, type="link")
}
mean.resid <- plyr::ddply(adults.count, "Year", summarize, mean.resid=mean(resid))
resid.fit <- lm( mean.resid ~ 1, data=mean.resid)
dwres1 <- car::durbinWatsonTest(resid.fit)
print(dwres1)
dwres2 <- lmtest::dwtest(resid.fit)
print(dwres2)
})
##
##
##
## *** Check for autocorrelation for species A-AMGR
## Warning in summary.lm(model): essentially perfect fit: summary may be
## unreliable
## lag Autocorrelation D-W Statistic p-value
## 1 -0.1666667 1.5 NA
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 1.5, p-value = 0.3328
## alternative hypothesis: true autocorrelation is greater than 0
##
##
##
##
## *** Check for autocorrelation for species A-AMMA
## Warning in summary.lm(model): essentially perfect fit: summary may be
## unreliable
## lag Autocorrelation D-W Statistic p-value
## 1 -0.1666667 1.5 NA
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 1.5, p-value = 0.3328
## alternative hypothesis: true autocorrelation is greater than 0
##
##
##
##
## *** Check for autocorrelation for species A-ANBO
## lag Autocorrelation D-W Statistic p-value
## 1 -0.5254895 2.576468 0.426
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 2.5765, p-value = 0.6954
## alternative hypothesis: true autocorrelation is greater than 0
##
##
##
##
## *** Check for autocorrelation for species A-PSRE
## lag Autocorrelation D-W Statistic p-value
## 1 -0.0006940065 1.002082 NA
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 1.0021, p-value = 0.01282
## alternative hypothesis: true autocorrelation is greater than 0
##
##
##
##
## *** Check for autocorrelation for species A-RAAU
## lag Autocorrelation D-W Statistic p-value
## 1 -0.0006940065 1.002082 0.06
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 1.0021, p-value = 0.01282
## alternative hypothesis: true autocorrelation is greater than 0
##
##
##
##
## *** Check for autocorrelation for species A-PRSE
## lag Autocorrelation D-W Statistic p-value
## 1 -0.6666667 3 NA
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 3, p-value = 1
## alternative hypothesis: true autocorrelation is greater than 0
##
##
##
##
## *** Check for autocorrelation for species ALL.species
## lag Autocorrelation D-W Statistic p-value
## 1 -0.02167572 1.065027 0.164
## Alternative hypothesis: rho != 0
##
## Durbin-Watson test
##
## data: resid.fit
## DW = 1.065, p-value = 0.1138
## alternative hypothesis: true autocorrelation is greater than 0
# extract a table of statistics for each study area x species combiations
adult.slopes <- ldply(fits, function(x){
adults.count <- x$fit$data
adult.fit <- x$fit
if(length(unique(adults.count$Transect.Label))==1){
adult.slopes <- data.frame(
Study.Area.Name =adults.count$Study.Area.Name[1],
slope = coef(adult.fit)["Year"],
slope.se = sqrt(diag(vcov(adult.fit)))["Year"],
p.value = summary(adult.fit)$coefficients[row.names(summary(adult.fit)$coefficients)=="Year" ,"Pr(>|t|)"],
stringsAsFactors=FALSE
)
}
if(length(unique(adults.count$Transect.Label))>1){
adult.slopes <- data.frame(
Study.Area.Name =adults.count$Study.Area.Name[1],
slope = fixef(adult.fit)["Year"],
slope.se = summary(adult.fit)$coefficients["Year","Pr(>|t|)"],
p.value = summary(adult.fit)$coefficients[row.names(summary(adult.fit)$coefficients)=="Year" ,"Pr(>|t|)"],
stringsAsFactors=FALSE)
}
adult.slopes
})
adult.slopes
## Species Study.Area.Name slope slope.se p.value
## 1 A-AMGR Alice Lake Park 0.000000e+00 1.224745e+00 1.000000e+00
## 2 A-AMMA Alice Lake Park 0.000000e+00 1.224745e+00 1.000000e+00
## 3 A-ANBO Alice Lake Park 1.513231e+00 5.285470e-01 2.139275e-01
## 4 A-PSRE Alice Lake Park 2.264197e+01 9.999992e-01 2.809852e-02
## 5 A-RAAU Alice Lake Park -2.264197e+01 9.999992e-01 2.809852e-02
## 6 A-PRSE Alice Lake Park 4.340178e-13 1.732051e+00 1.000000e+00
## 7 ALL.species Alice Lake Park 6.931472e-01 8.976133e-12 8.244113e-12
# compute the fitted values from the model
adult.fitted <- ldply(fits, function(x){
adults.count <- x$fit$data
adult.fit <- x$fit
Transect.Label <- unique(as.character(adults.count$Transect.Label))
newdata <- expand.grid(Year=seq(min(adults.count$Year, na.rm=TRUE),max(adults.count$Year, na.rm=TRUE), .1),
Transect.LabelF=Transect.Label)
newdata$Transect.Label <-as.character(newdata$Transect.LabelF)
if(length(unique(adults.count$Transect.Label))==1){
newdata$pred.mean <- predict(adult.fit, newdata=newdata,type="response")
}
if(length(unique(adults.count$Transect.Label))>1){
newdata$pred.mean <- exp(predict(adult.fit, newdata=newdata,type="response", re.form=~0))
}
# we now must average over all of the transect labels
adult.fitted <- plyr::ddply(newdata, c("Year","Transect.Label"), plyr::summarize, pred.mean=mean(pred.mean))
adult.fitted$Study.Area.Name <- adults.count$Study.Area.Name
adult.fitted
})
adult.fitted
## Species Year Transect.Label pred.mean Study.Area.Name
## 1 A-AMGR 2013.0 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 2 A-AMGR 2013.1 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 3 A-AMGR 2013.2 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 4 A-AMGR 2013.3 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 5 A-AMGR 2013.4 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 6 A-AMGR 2013.5 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 7 A-AMGR 2013.6 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 8 A-AMGR 2013.7 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 9 A-AMGR 2013.8 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 10 A-AMGR 2013.9 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 11 A-AMGR 2014.0 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 12 A-AMGR 2014.1 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 13 A-AMGR 2014.2 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 14 A-AMGR 2014.3 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 15 A-AMGR 2014.4 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 16 A-AMGR 2014.5 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 17 A-AMGR 2014.6 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 18 A-AMGR 2014.7 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 19 A-AMGR 2014.8 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 20 A-AMGR 2014.9 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 21 A-AMGR 2015.0 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 22 A-AMMA 2013.0 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 23 A-AMMA 2013.1 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 24 A-AMMA 2013.2 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 25 A-AMMA 2013.3 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 26 A-AMMA 2013.4 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 27 A-AMMA 2013.5 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 28 A-AMMA 2013.6 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 29 A-AMMA 2013.7 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 30 A-AMMA 2013.8 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 31 A-AMMA 2013.9 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 32 A-AMMA 2014.0 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 33 A-AMMA 2014.1 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 34 A-AMMA 2014.2 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 35 A-AMMA 2014.3 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 36 A-AMMA 2014.4 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 37 A-AMMA 2014.5 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 38 A-AMMA 2014.6 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 39 A-AMMA 2014.7 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 40 A-AMMA 2014.8 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 41 A-AMMA 2014.9 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 42 A-AMMA 2015.0 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 43 A-ANBO 2013.0 Fawn Lake Transect 1.528729e-01 Alice Lake Park
## 44 A-ANBO 2013.1 Fawn Lake Transect 1.778481e-01 Alice Lake Park
## 45 A-ANBO 2013.2 Fawn Lake Transect 2.069036e-01 Alice Lake Park
## 46 A-ANBO 2013.3 Fawn Lake Transect 2.407060e-01 Alice Lake Park
## 47 A-ANBO 2013.4 Fawn Lake Transect 2.800307e-01 Alice Lake Park
## 48 A-ANBO 2013.5 Fawn Lake Transect 3.257801e-01 Alice Lake Park
## 49 A-ANBO 2013.6 Fawn Lake Transect 3.790036e-01 Alice Lake Park
## 50 A-ANBO 2013.7 Fawn Lake Transect 4.409223e-01 Alice Lake Park
## 51 A-ANBO 2013.8 Fawn Lake Transect 5.129569e-01 Alice Lake Park
## 52 A-ANBO 2013.9 Fawn Lake Transect 5.967600e-01 Alice Lake Park
## 53 A-ANBO 2014.0 Fawn Lake Transect 6.942542e-01 Alice Lake Park
## 54 A-ANBO 2014.1 Fawn Lake Transect 8.076762e-01 Alice Lake Park
## 55 A-ANBO 2014.2 Fawn Lake Transect 9.396283e-01 Alice Lake Park
## 56 A-ANBO 2014.3 Fawn Lake Transect 1.093138e+00 Alice Lake Park
## 57 A-ANBO 2014.4 Fawn Lake Transect 1.271726e+00 Alice Lake Park
## 58 A-ANBO 2014.5 Fawn Lake Transect 1.479492e+00 Alice Lake Park
## 59 A-ANBO 2014.6 Fawn Lake Transect 1.721200e+00 Alice Lake Park
## 60 A-ANBO 2014.7 Fawn Lake Transect 2.002396e+00 Alice Lake Park
## 61 A-ANBO 2014.8 Fawn Lake Transect 2.329533e+00 Alice Lake Park
## 62 A-ANBO 2014.9 Fawn Lake Transect 2.710115e+00 Alice Lake Park
## 63 A-ANBO 2015.0 Fawn Lake Transect 3.152873e+00 Alice Lake Park
## 64 A-PSRE 2013.0 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 65 A-PSRE 2013.1 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 66 A-PSRE 2013.2 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 67 A-PSRE 2013.3 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 68 A-PSRE 2013.4 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 69 A-PSRE 2013.5 Fawn Lake Transect 1.778588e-15 Alice Lake Park
## 70 A-PSRE 2013.6 Fawn Lake Transect 1.711606e-14 Alice Lake Park
## 71 A-PSRE 2013.7 Fawn Lake Transect 1.647146e-13 Alice Lake Park
## 72 A-PSRE 2013.8 Fawn Lake Transect 1.585113e-12 Alice Lake Park
## 73 A-PSRE 2013.9 Fawn Lake Transect 1.525417e-11 Alice Lake Park
## 74 A-PSRE 2014.0 Fawn Lake Transect 1.467969e-10 Alice Lake Park
## 75 A-PSRE 2014.1 Fawn Lake Transect 1.412685e-09 Alice Lake Park
## 76 A-PSRE 2014.2 Fawn Lake Transect 1.359482e-08 Alice Lake Park
## 77 A-PSRE 2014.3 Fawn Lake Transect 1.308283e-07 Alice Lake Park
## 78 A-PSRE 2014.4 Fawn Lake Transect 1.259013e-06 Alice Lake Park
## 79 A-PSRE 2014.5 Fawn Lake Transect 1.211598e-05 Alice Lake Park
## 80 A-PSRE 2014.6 Fawn Lake Transect 1.165968e-04 Alice Lake Park
## 81 A-PSRE 2014.7 Fawn Lake Transect 1.122057e-03 Alice Lake Park
## 82 A-PSRE 2014.8 Fawn Lake Transect 1.079800e-02 Alice Lake Park
## 83 A-PSRE 2014.9 Fawn Lake Transect 1.039134e-01 Alice Lake Park
## 84 A-PSRE 2015.0 Fawn Lake Transect 1.000000e+00 Alice Lake Park
## 85 A-RAAU 2013.0 Fawn Lake Transect 1.000000e+00 Alice Lake Park
## 86 A-RAAU 2013.1 Fawn Lake Transect 1.039134e-01 Alice Lake Park
## 87 A-RAAU 2013.2 Fawn Lake Transect 1.079800e-02 Alice Lake Park
## 88 A-RAAU 2013.3 Fawn Lake Transect 1.122057e-03 Alice Lake Park
## 89 A-RAAU 2013.4 Fawn Lake Transect 1.165968e-04 Alice Lake Park
## 90 A-RAAU 2013.5 Fawn Lake Transect 1.211598e-05 Alice Lake Park
## 91 A-RAAU 2013.6 Fawn Lake Transect 1.259013e-06 Alice Lake Park
## 92 A-RAAU 2013.7 Fawn Lake Transect 1.308283e-07 Alice Lake Park
## 93 A-RAAU 2013.8 Fawn Lake Transect 1.359482e-08 Alice Lake Park
## 94 A-RAAU 2013.9 Fawn Lake Transect 1.412685e-09 Alice Lake Park
## 95 A-RAAU 2014.0 Fawn Lake Transect 1.467969e-10 Alice Lake Park
## 96 A-RAAU 2014.1 Fawn Lake Transect 1.525417e-11 Alice Lake Park
## 97 A-RAAU 2014.2 Fawn Lake Transect 1.585113e-12 Alice Lake Park
## 98 A-RAAU 2014.3 Fawn Lake Transect 1.647146e-13 Alice Lake Park
## 99 A-RAAU 2014.4 Fawn Lake Transect 1.711606e-14 Alice Lake Park
## 100 A-RAAU 2014.5 Fawn Lake Transect 1.778588e-15 Alice Lake Park
## 101 A-RAAU 2014.6 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 102 A-RAAU 2014.7 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 103 A-RAAU 2014.8 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 104 A-RAAU 2014.9 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 105 A-RAAU 2015.0 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 106 A-PRSE 2013.0 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 107 A-PRSE 2013.1 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 108 A-PRSE 2013.2 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 109 A-PRSE 2013.3 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 110 A-PRSE 2013.4 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 111 A-PRSE 2013.5 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 112 A-PRSE 2013.6 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 113 A-PRSE 2013.7 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 114 A-PRSE 2013.8 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 115 A-PRSE 2013.9 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 116 A-PRSE 2014.0 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 117 A-PRSE 2014.1 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 118 A-PRSE 2014.2 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 119 A-PRSE 2014.3 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 120 A-PRSE 2014.4 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 121 A-PRSE 2014.5 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 122 A-PRSE 2014.6 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 123 A-PRSE 2014.7 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 124 A-PRSE 2014.8 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 125 A-PRSE 2014.9 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 126 A-PRSE 2015.0 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 127 ALL.species 2013.0 Fawn Lake Transect 1.000000e+00 Alice Lake Park
## 128 ALL.species 2013.1 Fawn Lake Transect 1.071773e+00 Alice Lake Park
## 129 ALL.species 2013.2 Fawn Lake Transect 1.148698e+00 Alice Lake Park
## 130 ALL.species 2013.3 Fawn Lake Transect 1.231144e+00 Alice Lake Park
## 131 ALL.species 2013.4 Fawn Lake Transect 1.319508e+00 Alice Lake Park
## 132 ALL.species 2013.5 Fawn Lake Transect 1.414214e+00 Alice Lake Park
## 133 ALL.species 2013.6 Fawn Lake Transect 1.515717e+00 Alice Lake Park
## 134 ALL.species 2013.7 Fawn Lake Transect 1.624505e+00 Alice Lake Park
## 135 ALL.species 2013.8 Fawn Lake Transect 1.741101e+00 Alice Lake Park
## 136 ALL.species 2013.9 Fawn Lake Transect 1.866066e+00 Alice Lake Park
## 137 ALL.species 2014.0 Fawn Lake Transect 2.000000e+00 Alice Lake Park
## 138 ALL.species 2014.1 Fawn Lake Transect 2.143547e+00 Alice Lake Park
## 139 ALL.species 2014.2 Fawn Lake Transect 2.297397e+00 Alice Lake Park
## 140 ALL.species 2014.3 Fawn Lake Transect 2.462289e+00 Alice Lake Park
## 141 ALL.species 2014.4 Fawn Lake Transect 2.639016e+00 Alice Lake Park
## 142 ALL.species 2014.5 Fawn Lake Transect 2.828427e+00 Alice Lake Park
## 143 ALL.species 2014.6 Fawn Lake Transect 3.031433e+00 Alice Lake Park
## 144 ALL.species 2014.7 Fawn Lake Transect 3.249010e+00 Alice Lake Park
## 145 ALL.species 2014.8 Fawn Lake Transect 3.482202e+00 Alice Lake Park
## 146 ALL.species 2014.9 Fawn Lake Transect 3.732132e+00 Alice Lake Park
## 147 ALL.species 2015.0 Fawn Lake Transect 4.000000e+00 Alice Lake Park
# Plot with trend line
adult.plot.summary <- ggplot2::ggplot(data=adults.count,
aes(x=Year, y=n.adults))+
ggtitle("Amphibian adults count ")+
ylab("Amphibian adults Count")+
geom_point(size=3, aes(color=Transect.Label))+
facet_wrap(~interaction(Study.Area.Name,Species), ncol=2, scales="free_y")+
geom_line(data=adult.fitted, aes(x=Year,y=pred.mean))+
scale_x_continuous(breaks=min(adults.count$Year,na.rm=TRUE):max(adults.count$Year,na.rm=TRUE))+
geom_text(data=adult.slopes, aes(x=min(adults.count$Year, na.rm=TRUE), y=max(adults.count$n.adults, na.rm=TRUE)),
label=paste("Slope (on log scale) : ",round(adult.slopes$slope,2),
" ( SE " ,round(adult.slopes$slope.se,2),")",
" \np :" ,round(adult.slopes$p.value,3)),
hjust="left",vjust="top")+
theme(legend.position="top")
adult.plot.summary
ggsave(plot=adult.plot.summary,
file=paste(file.prefix,'-adult-plot-summary.png',sep=""),
h=6, w=6, units="in", dpi=300)